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Executive  Summary 


The  Segmented  Mirror  Telescope  (SMT)  was  built  to  develop  new  techniques  for  rapid  con¬ 
struction  of  optical  space-based  telescopes.  This  system  has  been  transferred  to  the  Naval 
Postgraduate  School  (NPS)  and  is  now  hosted  by  the  Spacecraft  Research  and  Design  Cen¬ 
ter  (SRDC).  As  a  part  of  this  initiative,  the  SRDC  now  is  the  Adaptive  Optics  Center  of  Ex¬ 
cellence  for  National  Security  and  is  actively  pursuing  research  focused  on  the  SMT  and  other 
optics  systems. 

The  most  complex  problem  of  optical  systems  is  constructing  large  primary  mirrors  of  sufficient 
optical  quality.  Traditional  monolithic  mirrors  can  take  up  to  four  years  to  construct  and  test.  If 
the  mirror  is  segmented  and  split  into  several  smaller  mirrors,  the  process  becomes  much  more 
rapid.  The  mirror  segments  for  the  SMT  were  constructed  in  approximately  18  months.  These 
smaller  mirrors  can  also  be  combined  to  form  much  larger  mirrors  than  could  have  been  built 
by  the  traditional  monolithic  processes. 

The  trade-off  in  this  new  technique  is  an  increased  complexity  of  the  system.  The  monolithic 
mirror  is  built  to  act  as  a  coherent  reflector  across  the  entire  mirror  surface.  The  larger  size 
also  means  that  it  is  not  as  susceptible  to  damage  during  the  rocket  launch  into  orbit.  For  the 
segmented  system,  the  mirrors  must  be  aligned  correctly.  If  this  process  is  not  done  properly  and 
the  mirror  height  between  segments  varies  too  much,  the  resulting  imaging  is  worse.  Phasing 
of  the  telescope,  the  highly  sensitive  alignment  process,  must  be  done  to  ensure  this  does  not 
occur. 

In  addition  to  the  mirror  surface  being  segmented,  the  rear  surface  also  has  surface-parallel 
piezoelectric  actuators.  These  actuators  can  correct  for  mirror  imperfections.  This  arrangement 
is  referred  as  an  active-optics  deformable  mirror  system.  This  allows  for  increased  variance 
tolerances  in  the  mirror  manufacturing  process  since  the  mirror  can  be  corrected  when  the  ac¬ 
tuators  are  energized.  This  technique  lowers  costs  associated  with  construction  and  testing  but 
adds  complexity  to  the  control  system. 

In  contrast  to  the  active-optics  system,  SMT  is  an  adaptive-optics  system.  The  telescope  in¬ 
cludes  a  Shack-Hartmann  Wavefront  Sensor  (WFS),  which  provides  information  about  the  qual¬ 
ity  of  the  light  reaching  the  camera.  The  light  may  lose  phase  coherency  and  become  less  planar. 
This  occurs  from  either  atmospheric  disturbances,  structural  jitter  and/or  mirror  imperfections. 
Each  source  of  error  is  summed  to  provide  the  total  error.  The  WFS  provides  the  feedback 


mechanism  for  the  mirror  actuators.  A  closed  control  loop  would  be  suited  for  this  situation 
to  minimize  wavefront  error  by  correcting  the  light  to  have  a  planar  wavefront  at  the  camera 
detector. 

In  order  to  correct  the  light,  an  estimate  of  the  wavefront  must  be  formed.  An  implementation 
of  the  Poyneer  algorithm  is  developed  in  MATLAB  to  handle  other  aperture  shapes,  such  as 
the  hexagons  used  in  the  SMT.  The  sampling  geometry  for  the  Shack-Hartmann  is  applied 
to  the  algorithm.  The  algorithm  provides  the  information  necessary  to  be  used  in  a  real-time 
implementation  of  the  wavefront  reconstruction  such  that  it  could  be  used  in  a  feedback  control 
loop. 

A  novel  controller  design  that  uses  the  two-dimensional  fast  Fourier  transform  (FFT2)  is  pro¬ 
posed  in  this  thesis  to  solve  this  problem.  First,  a  simple  model  is  developed  where  a  deformable 
mirror  surface  is  treated  as  a  grid  of  masses  and  springs.  The  actuators  have  surface  normal 
forces  applied  to  individual  masses.  The  equation  to  describe  this  arrangement  is  a  triple  convo¬ 
lution  in  both  the  space  and  time  domains.  Because  of  the  complexity  of  the  coupled  equations, 
the  FFT2  is  employed.  In  the  frequency  domain,  the  problem  transforms  into  a  much  simpler 
arrangement  of  many  uncoupled,  independent  equations. 

After  developing  the  mirror  model,  we  designed  the  controller  model.  Instead  of  the  traditional 
design  of  one  controller  for  each  actuator,  one  controller  for  each  spatial  mode  of  actuators 
is  used.  Each  controller  then  provides  actuator  commands  for  every  actuator.  The  inverse 
transform  is  performed  to  recover  the  specific  commands  for  each  actuator. 

The  mirror  model  and  the  controller  are  implemented  in  this  thesis  using  MATLAB  Simulink. 
The  simulation  takes  advantage  of  the  powerful  discrete  filter  blocksets  to  efficiently  process 
the  transfer  functions  in  parallel.  The  mirror  aperture  size  is  set  to  be  256x256  to  show  the 
algorithm  is  scalable  to  much  larger  apertures  than  are  currently  used  by  state-of-the-art  optical 
systems. 

A  closed-loop  feedback  control  simulation  is  performed  with  time-varying  phase  wavefront 
information.  The  simulation  is  run  to  verify  the  design  of  the  controller  and  the  wavefront  error 
decrease  as  the  actuators  affect  the  mirror  surface.  The  results  are  presented  in  the  thesis.  This 
form  of  a  controller  can  be  used  on  the  actual  SMT  after  determining  the  actual  coefficients  of 
the  mechanical  structure  of  the  system. 

The  implementation  of  a  discrete  Fourier  transform  (DFT)  wavefront  reconstruction  and  de- 


xiv 


formable  mirror  model  and  controller  is  developed  in  this  thesis.  These  tools  will  be  useful 
in  the  future  research  of  the  SMT  and  can  be  applied  to  an  actual  system  for  experimental  re¬ 
sults.  Any  techniques  that  improve  the  imaging  quality  of  the  SMT  can  be  applied  to  other 
telescope  systems  to  ensure  this  technology  continues  to  progress  with  the  goal  of  providing 
higher  quality  image  products  for  less  cost. 
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CHAPTER  1: 

Introduction 


Traditional  space-based  optical  telescopes  employ  a  single,  large  primary  mirror  structure.  This 
mirror  has  to  be  painstakingly  manufactured  to  minimize  deformation  errors  across  the  surface. 
The  testing  may  take  months  to  years,  depending  on  the  program  requirements.  These  large 
mirrors  are  costly  to  manufacture,  test  and  schedule  for  spacecraft.  The  astronomy  community 
has  developed  a  number  of  technologies  that  allow  for  segmented  mirrors  to  replace  single,  large 
mirrors.  This  technology  is  now  being  developed  by  the  space-based  telescope  community  for 
the  next  generation  of  telescopes  [1]. 

The  Naval  Postgraduate  School  (NPS)  hosts  a  Segmented  Mirror  Telescope  (SMT),  originally 
produced  by  the  National  Reconnaissance  Office  (NRO)  as  a  proof-of-concept  experimental 
spacecraft  telescope  [2].  It  incorporates  many  aggressive  concepts  into  a  single  working  unit  to 
test  and  leam  about  the  manufacturing  of  such  an  instrument. 

The  primary  mirror  is  segmented  into  six  equal  hexagons  arranged  in  the  same  manner  as  a 
traditional  mirror.  The  segments  are  paired  with  each  other  and  have  two  sets  of  hinges.  In  this 
arrangement,  the  segmented  mirror  can  be  folded  up  for  launch  and  mechanically  unfold  and 
effloresce  the  segments.  This  allows  for  a  smaller  launch  configuration,  which  increases  options 
for  launch  vehicles. 

In  addition  to  the  primary  mirror  being  segmented,  it  also  has  actuators  across  the  rear  of  the 


Figure  1.1:  SMT  segmented  mirror  orientation. 
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surface.  There  are  three  types  of  actuators:  coarse,  fine  and  face  sheet.  These  actuators  can 
modify  the  mirror  surface  in  such  a  way  as  to  correct  for  abnormalities  in  the  approaching  light. 
The  coarse  and  fine  actuators  are  used  as  a  calibration  to  give  the  primary  mirror  a  high-quality 
optical  surface.  These  actuators  make  the  system  an  active-optics  system.  With  a  wavefront 
sensor  these  distortions  are  measured,  and  a  control  system  adjusts  the  face-sheet  actuators  to 
compensate.  In  this  manner,  the  telescope  can  also  be  considered  an  adaptive -optics  system, 
which  is  an  improvement  over  active  optics.  What  makes  this  unique  is  that  the  primary  mirror 
has  actuators,  whereas  most  adaptive-optics  systems  make  use  of  a  smaller  deformable  mirror 
such  as  a  fold  mirror  that  is  after  the  primary  in  the  optical  path  [3]. 

The  National  Security  Space  Strategy  [4]  outlines  the  importance  of  systems  like  the  SMT. 
Similar  systems  offer  “unprecedented  advantages  in  national  decision-making,  military  oper¬ 
ations  and  homeland  security.”  The  United  States  must  maintain  the  benefits  afforded  by  our 
exquisite  systems  in  the  evolving  strategic  environment.  Developing  technology  that  improves 
the  design  of  the  SMT  directly  supports  these  goals. 


1.1  Purpose 

Some  of  the  background  fundamental  concepts  required  for  a  space-based  optical  telescope  that 
functions  as  an  adaptive-optics  system  are  developed  in  this  thesis.  The  SMT  has  a  unique 
geometry  that  requires  finesse  in  application  of  these  concepts.  The  concepts  are  introduced 
for  a  simple  square  mirror  and  then  applied  to  the  SMT  and  modeled.  Simulation  results  are 
presented.  The  improvements  made  to  the  SMT  model  and  controller  will  result  in  more  accu¬ 
rate  analysis  in  future  space  telescope  research.  This  work  benefits  government  agencies  and 
companies  interested  in  space  telescopes. 

The  purpose  of  this  research  is  to  improve  the  quality  of  the  imaging  results  of  the  telescope.  The 
ideal  wavefront  is  planar  at  the  camera  sensor;  although,  the  wavefront  can  be  distorted  from 
the  ideal  geometry  for  a  number  of  reasons.  Atmospheric  mixing  and  motion  can  cause  non- 
uniform  delays  in  the  wavefront,  a  common  problem  for  ground-based  astronomical  telescopes. 
For  spacecraft  telescopes,  the  structural  jitter  from  the  motion  of  changing  the  satellite  pointing 
direction  can  deform  the  optics.  For  Earth-observation  satellites,  this  is  a  great  concern  because 
the  pointing  direction  is  changed  often  during  operation.  In  this  research,  a  sample  wavefront 
is  used  in  the  discussion.  The  three-dimensional  image  of  the  wavefront  is  shown  in  Figure  1.2. 
The  red  areas  indicate  the  leading  wavefront,  and  the  blue  areas  indicate  the  lagging  wavefront. 
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Figure  1.2:  Projection  of  the  sample  wavefront. 


1.2  Overview 

The  background  theory  of  wavefront  reconstruction,  which  is  a  process  employed  to  correct  the 
optical  system,  is  provided  in  Chapter  2.  This  process  is  applied  to  a  square  and  circular  aper¬ 
tures,  which  was  developed  in  previous  work  [5,6].  Wavefront  reconstruction  is  then  discussed 
for  the  hexagonal  aperture  and  the  unique  characteristics  of  the  SMT. 

The  two-dimensional  fast  Fourier  transform  (FFT2)  is  used  to  model  the  SMT  in  Chapter  3. 
Use  of  the  FFT2  leads  to  a  computationally  efficient  algorithm  that  yields  good  results.  This 
technique  yields  decoupled  dynamics  and  a  more  efficient  control  design.  This  improves  the 
scalability  of  the  processing  such  that  larger  apertures  can  be  computed  without  large  increases 
in  computational  time. 

The  implementation  and  simulation  work  done  for  this  thesis  is  covered  in  Chapter  4.  Wavefront 
reconstruction  from  example  generated  data  is  presented.  The  technique  can  also  be  applied  to 
noisy  slope  data  from  actual  systems.  The  Deformable  Mirror  (DM)  controller  using  the  discrete 
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Fourier  transform  (DFT)  is  also  described,  and  the  technique  of  building  the  actual  model  and 
its  simulation  are  also  presented.  The  transfer  function  equivalents  are  considered  to  reduce  the 
order  of  the  system.  Sample  phase  data  is  used  to  illustrate  the  performance  of  the  system. 

The  summary  of  the  work  and  conclusions  with  recommended  follow  on  research  for  the  SMT 
are  discussed  in  Chapter  5. 

1.3  Research  Contributions 

Square  and  circular  apertures  are  discussed  in  existing  literature  and  analyzed  in  great  detail 
[5-8].  Although  upcoming  telescopes  use  hexagonal  apertures,  they  are  not  widely  treated 
in  the  literature.  The  required  modifications  for  the  wavefront  reconstruction  algorithm  of  a 
circular  aperture  to  correctly  model  a  hexagonal  aperture  are  provided. 

Zernike  polynomials  are  used  commonly  in  literature  and  practice  to  characterize  systems  and 
wavefronts  [9].  Although  the  polynomials  are  particularly  suitable  to  model  optical  wavefronts, 
they  are  computationally  inefficient.  An  alternate  technique  that  is  computationally  efficient 
is  desired  for  realtime,  closed-loop  operations  of  a  DM.  A  spatial  frequency  wavefront  recon¬ 
struction  technique  is  developed  as  a  viable  alternate  to  Zemike  polynomials. 

A  further  contribution  is  also  made  in  the  control  system  design.  In  a  traditional  control  design, 
commands  are  generated  for  each  actuator  individually  in  the  system.  This  does  not  scale 
efficiently  as  the  number  of  actuators  increases.  A  modal  controller  is  developed  to  correct  for 
wavefront  distortions  in  the  spatial  frequency  domain. 


4 


CHAPTER  2: 

Wavefront  Reconstruction 


2.1  Theory 

Wavefront  reconstruction  is  the  mathematical  process  of  using  the  sensed  wavefront  gradient 
information  and  determining  an  estimated  equivalent  wavefront.  An  approaching  wavefront  has 
very  fine  perturbations  in  the  phase  which  can  be  observed  in  the  resulting  images  as  distortions. 
The  differences  in  phase  can  be  specified  by  using  the  optical  path  length  [10] 

2k 

8(f)  —  —n8d  (2.1) 

A 

where  8(j)  is  the  difference  in  phase,  A  is  the  wavelength  of  light,  n  is  the  index  of  refraction 
of  the  media  and  8d  is  the  minute  difference  in  the  physical  path  length,  expressed  as  a  lump 
sum  of  all  of  the  distortions.  The  quantity  nd  is  the  optical  path  length.  In  the  laboratory 
environment,  the  medium  is  air.  The  optical  path  length  is  usually  considered  for  the  center 
wavelength  of  interest,  in  this  case  within  the  visible  light  band.  Thus,  the  correction  cannot  be 
applied  uniformly  for  all  wavelengths,  as  each  has  a  different  optical  path  length. 

These  differences  reduce  the  resolution  and  clarity  on  the  resulting  image.  To  correct  for  the 
phase  differences,  the  mirror  is  properly  reformed  by  a  set  of  actuators.  This  can  be  seen  in 
Figure  1.2,  where  the  mirror  needs  to  provide  a  negative  phase  adjustment  for  the  red  areas  by 
pushing  the  mirror  back  and  a  positive  phase  adjustment  for  the  blue  areas  by  pushing  the  mirror 
forward,  so  that  when  the  light  arrives  at  the  mirror  surface,  the  wavefront  reflection  is  planar. 
The  end  result  is  improved  imaging. 

In  the  actual  implementation  of  Equation  (2.1),  the  phase  (j)  is  an  estimate  that  cannot  fully 
eliminate  all  errors.  The  phase  estimates  (j)(x,y)  are  computed  across  a  two-dimensional  grid  of 
values  on  the  basis  of  observed  local  gradients.  The  gradients  are  sampled  at  a  low  spatial  fre¬ 
quency  across  the  aperture  and,  as  a  consequence  of  the  sampling,  can  only  reconstruct  specific 
low-order  modes  of  the  wavefront.  These  gradients  are  then  processed  to  determine  a  wave- 
front  that  has  local  gradients  of  the  same  form  as  the  original  wave.  The  device  that  samples 
the  wavefront  is  called  the  wavefront  sensor  [11].  There  are  several  different  wavefront  sensors; 
the  one  that  is  discussed  in  this  thesis  is  the  Shack-Hartmann  sensor  [12],  which  is  a  grid  of 
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Figure  2.1:  Square  grid  of  phase  points  that  contribute  to  gradients  in  (a)  H udgi n  geometry 
and  (b)  Fried  geometry.  Shack-Hartmann  sensors  are  shown  in  red  dashed  lines  with  hexagon 
lenslets. 


small  lenses  that  are  side-by-side.  Each  lenslet  focuses  the  received  light  to  a  point  on  the  focal 
plane  of  a  secondary  camera.  If  the  point  does  not  appear  at  the  center  of  the  lenslet  optical  axis 
on  the  camera  focal  plane,  then  the  light  has  wavefront  variation.  The  position  of  these  points 
relative  to  their  expected  centered  position  yields  the  information  on  the  phase  gradient. 

There  are  two  common  sensor  geometries  developed  by  Hudgin  [7]  and  Fried  [8]  shown  in 
Figure  2.1.  These  geometries  vary  by  the  location  of  the  Wavefront  Sensor  (WFS)  relative 
to  the  phase  points,  which  are  abstract  locations  in  the  optical  plane  that  have  the  equivalent 
gradients  between  adjacent  phase  points  as  registered  by  the  sensors.  The  Hudgin  geometry 
is  easier  to  work  with  initially  in  simulation  but  has  the  drawback  that  the  sensor’s  lenslets 
have  overlap,  which  is  not  commonly  done  in  actual  implementation.  Fried  geometry  is  more 
commonly  used  to  model  the  actual  Shack-Hartmann  sensors.  The  SMT  uses  hexagon-shaped 
lenslets  for  its  sensors.  These  pack  together  densely  without  overlap  and  focus  all  available 
light  on  the  sensor  detection  focal  plane. 

Wavefronts  can  be  reconstructed  using  a  variety  of  basis  functions.  Some  reconstructions  use 
Zernike  polynomials,  which  are  particularly  suitable  to  describe  the  common  optical  character¬ 
izations  of  astigmatism,  coma,  defocus  and  others  [9].  Although  a  full  expansion  is  computa¬ 
tionally  inefficient,  only  a  few  terms  are  needed  to  characterize  a  wavefront. 

An  alternate  decomposition  by  the  DFT  yields  the  desirable  property  of  being  computationally 
efficient  and  still  only  providing  a  few  dominant  coefficients  [6].  The  work  in  this  thesis  uses 
this  technique. 
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2.2  Square  Aperture 


The  square  aperture  is  a  primitive  geometry  that  allows  for  a  simplified  reconstruction  algo¬ 
rithm.  The  design  assumes  that  the  actuators  and  phase  points  are  co-located.  Thus,  the  actu¬ 
ators  directly  influence  the  phase  points  (and,  thereby,  the  gradients).  For  a  NxN  grid  of  phase 
points  (f>  \m.  n\.  Hudgin’s  geometry  model  defines  the  gradients  as  the  difference  of  the  neighbor 
phase  points: 


sx  [m,  n]  =  (j)  [m  +  1 ,  n]  —  (j)  [m,  n]  (2.2a) 

Sy  [m,  n\  =  (j)  [m,  n  +  1]  —  (j)  [m,  n ] .  (2.2b) 

Alternatively,  Fried’s  model  defines  the  gradients  as  the  average  of  the  neighbor  phase  points: 

sx[m,n]  —  ~((j)[m+  1  ,n]  —  (j)[m.n]  +<j)[m+  1  ,n+  1]  —  (/)[m,n+  1])  (2.3a) 

sy[m,n]  —  ~((j)[m,n  +  1]  —  (j)[m,n\  +  0[m+  l,n  +  1]  —  (j)[m+  l,n]).  (2.3b) 

In  either  case,  the  result  is  two  sets  of  slope  data  that  are  an  ( N  —  l)x(N  —  1)  grid  of  values, 
rather  than  NxN  because  the  slopes  can  only  be  computed  between  data  points.  Slopes  are 
defined  as  the  local  gradients  separated  into  the  individual  components,  in  this  case  sx  and  sy. 
Using  the  notion  that  any  closed  path  along  the  mirror  surface  must  result  in  a  net  slope  change 
of  zero,  Freischlad  determined  the  last  index  of  values  along  the  boundaries  as  [5]: 

N- 1 

sx[N,n\  =  -  Y,  sx[k,n\  (2.4a) 

k=\ 

N- 1 

sy[m,N]  =  -  ^  sy[m,k].  (2.4b) 

k=  1 


The  two  sets  of  slope  data  are  now  NxN.  An  alternate  method  is  to  assume  that  (j)\in.n]  is 
periodic  and  that  (j)  [1,  N+ 1]  =  (j)  [1, 1]  and  so  on.  In  this  method,  the  last  slopes  can  be  computed, 
also  resulting  in  NxN  data.  This  periodic  assumption  is  included  in  the  definition  of  the  DFT. 

After  sx  and  sy  are  determined  from  the  sensors,  the  reconstruction  of  the  0  can  be  completed 
using  the  FFT2.  The  first  step  is  to  take  the  FFT2  of  the  sx  and  sy  data,  as 
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N-lN-l 

S*[M=  E  Y, 

(2.5a) 

p= 0  q= 0 

TV-t/V-l 

SyM=  E  E  hM^{kq+,q) ■ 

(2.5b) 

p= 0  q= 0 


The  data  can  be  filtered  in  the  frequency  domain.  This  filtering  is  equivalent  to  solving  Equa¬ 
tions  (2.3a)  and  (2.3b)  for  (j)  in  the  frequency  domain  to  obtain 


4>M 


o, 

j2nk  j2nl 

[e~-N--l]Sx{kJ]+{e—N-~l}Sy{k,l] 

[4  (sin2(f)+sin2{§))\ 


ifk,l=0 

otherwise 


j2xk 


j2nl  j2nl 

■l][e-^r+l]Sx[k,l]+[e—n-- 


j2nk 


[8  [sin2  ( ^  )cos2  ( f )  +sin2  ( f )  cos2  ( ^ ) )] 


if  k,l  —  0',k,l  —N/2 
otherwise. 


(2.6a) 


(2.6b) 


Equation  (2.6a)  is  valid  for  Hudgin  geometry,  while  Equation  (2.6b)  is  valid  for  Fried  geometry. 
The  resulting  wavefront  phase  estimate  is  available  after  taking  the  inverse  FFT2: 

#K»]  =  4  E'  e' ' (2.7) 

p= o  q= 0 

The  phase  value  at  the  phase  points  is  now  known.  This  information  can  then  be  processed  by 
the  controller  to  properly  deform  the  mirror  surface  and  improve  the  resulting  image.  Since  the 
traditional  deformable  mirror  has  the  actuators  co-located  with  the  phase  points,  the  controller 
can  easily  make  the  necessary  changes  that  directly  improves  the  imaging  results. 

2.2.1  Boundary  Conditions 

The  combination  of  the  sums  in  Equations  (2.4a)  and  (2.4b)  imply  that  there  are  artifacts  along 
the  boundary  of  two  of  the  four  sides.  The  sums  of  noisy  data  points  result  in  a  larger  variance 
of  the  noise  on  the  sum.  This  is  an  undesirable  consequence  of  this  technique. 
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The  boundary  is  where  the  majority  of  the  error  occurs  when  computing  the  wavefront.  These 
errors  can  result  in  large  discontinuities  which  affect  the  reconstruction  algorithm  at  other  points 
across  the  entire  aperture.  This  can  be  explained  since  the  impulse  response  of  the  algorithm  is 
shift  variant,  although  the  reconstruction  is  unbiased  [5].  That  is  to  say,  the  sums  in  Equations 
(2.4a)  and  (2.4b)  have  detrimental  effects  across  the  entire  aperture  in  reconstruction;  a  better 
technique  must  address  this  issue. 

2.3  Circular  Aperture 

The  circular  aperture  was  analyzed  by  Poyneer  [6] .  Wavefront  reconstruction  in  the  center  of 
the  aperture  is  relatively  straightforward  from  the  square  aperture  discussion  in  Section  2.2. 
However,  along  the  circular  boundary,  there  are  large  errors  in  the  slopes  as  calculated  by  Equa¬ 
tion  (2.3).  These  occur  because  the  area  external  to  the  circular  aperture  is  set  to  zero  in  the 
matrix  grid  of  data.  The  slope  data  is  then  calculated  incorrectly,  and  the  error  propagates  into 
the  reconstructed  wavefront.  Several  modifications  are  needed  to  understand  how  the  slopes  are 
affected  by  the  boundary. 

In  Figure  2.2,  the  gradients  are  shown  as  scaled  vectors  sampled  at  a  very  low  rate.  The  vectors 
near  the  boundary  contain  large  errors,  as  can  be  seen  by  the  larger  magnitude  and  even  opposite 
phase.  Additionally,  some  of  the  gradients  are  non-zero  outside  of  the  aperture. 

2.3.1  Boundary  Conditions 

The  definition  of  the  slopes  in  Equations  (2.2)  has  to  be  altered  to  accommodate  the  circular 
boundary.  At  each  boundary  point  [m,n],  let 

sx[m,  n]  =  ix[m.  n\  +  bx[m,  n]  (2.8a) 

sy  [m,  n\  =  iy  [, m ,  n\  +  by  [m,  n\  (2.8b) 

where  the  values  ix[m,n]  and  iy[m.n]  are  internal  slopes  near  the  boundary  that  are  measured 
using  Equations  (2.2);  these  are  the  necessary  values  to  compute  wave-front  phase  reconstruc¬ 
tion  correctly.  Additionally,  the  bx[m.n]  and  by [m.n]  are  the  boundary  slopes  that  are  incorrect; 
removal  of  these  terms  is  necessary  for  proper  wavefront  reconstruction. 

The  relationship  between  the  internal  and  boundary  gradients  are  defined  by  the  matrix 


M«  =  c 


(2.9) 
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Figure  2.2:  Overhead  view  of  gradients  for  sample  wavefront. 


where  M  is  a  matrix  determined  by  the  geometry,  u  is  a  vector  of  unknown  boundary  slopes  and 
c  is  a  vector  of  the  measured  slopes  which  slopes  cross  the  boundary. 

The  matrix  M  is  developed  from  the  equations  that  relate  the  entries  of  u  to  c.  To  solve  for  u, 
the  inverse  of  M  is  required.  Equation  (2.9)  has  an  infinite  number  of  solutions  since  it  does 
not  solve  for  the  DC  term  (colloquially  known  as  the  piston).  This  occurs  because,  while  the 
slopes  of  the  phase  are  known,  the  actual  phase  is  not.  Since  M  is  singular,  M  might  be  better 
conditioned  by  singular  value  decomposition  while  keeping  the  dominant  singular  values. 

In  Figure  2.3,  the  unknown  boundary  slopes  are  labeled  sequentially  in  u,  while  the  known 
internal  slopes  are  sx[m.  n ]  and  sy  [m.  n] .  The  u  slopes  all  cross  the  circular  boundary  edge  and  is 
shown  as  a  dotted  line.  The  zero  vectors  shown  indicate  that  the  two  connected  points  are  both 
external  to  the  circular  aperture. 

By  writing  the  mesh  equations  for  small  square  edges  between  adjacent  phase  points,  we  obtain 
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Figure  2.3:  Circular  aperture  showing  slopes  near  the  boundary. 


a  set  of  linear  equations.  This  method  generates  the  linear  combinations  necessary  to  solve  for 
the  non-noisy  data  exact  solution  or  the  noisy  data  least-squares  solution.  These  equations  di¬ 
rectly  relate  the  measured  slopes  to  the  unknowns  and  zero  external  slopes.  A  starting  location 
must  be  selected  along  the  boundary,  and  then  a  path  is  traced  along  the  boundary  until  the 
aperture  edge  is  closed.  Depending  on  the  geometry,  we  see  that  the  c  entry  is  a  linear  combi¬ 
nation  of  some  sx[m,n ]  and  sy[m,n].  Upon  solving  for  the  u  solution,  we  obtain  each  resulting 
u  entry  as  either  a  bx[m.n\  or  by\m,n\.  Knowing  the  bx[m,n\  or  by[m,n]  information,  we  solve 
Equation  (2.8)  for  the  needed  ix[m,n\  and  iy[m,n],  which  is  trivial.  The  results  correct  the  slope 
data  along  the  boundary,  and  the  data  is  now  used  for  reconstruction  of  the  phase  0. 
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Figure  2.4:  SMT  wavefront  sensor  topology  for  a  single  segment  with  modified  Fried  geometry 
phase  point  locations. 

2.4  Hexagonal  Aperture 

The  hexagonal  aperture  is  a  modification  of  a  circular  aperture.  It  is  more  desirable  for  multiple 
mirror  systems  than  the  latter  since  its  segments  can  easily  fit  together  in  a  continuous  fashion. 
It  has  been  commonly  used  on  many  recent  space  telescopes,  including  the  SMT  and  the  James 
Webb  Space  Telescope  (JWST). 

The  linear  edges  make  the  mirror  segments  easier  to  align  with  one  another,  such  as  by  using 
edge-detection  sensors.  These  sensors  are  able  to  detect  the  neighbor  segment  mirror  heights 
and  adjust  the  coarse  and  fine  actuators  to  place  the  heights  of  the  mirror  surfaces  very  closely  to 
one  another  at  the  nanometer  scale.  This  calibration  is  referred  to  as  phasing  the  telescope  and  is 
seldom  needed  to  be  repeated  for  on-orbit  operation  but  is  performed  often  for  the  ground-based 
telescopes. 

For  the  SMT,  the  61  WFS  per  segment  pictured  in  Figure  2.4  are  packed  together  densely. 
However,  their  combined  edge  does  not  result  in  an  exact  hexagonal  shape  along  the  outer  edge. 
The  phase  points  should  be  placed  along  the  edges  of  the  individual  lenslets.  This  makes  the 
grid  uniform  in  rectangular  coordinates.  This  is  a  more  convenient  arrangement  than  the  one 
originally  depicted  in  Figure  2.1.  If  the  phase  points  were  completely  inside  the  lenslet  for 
each  lenslet,  there  would  be  no  information  about  the  relationship  between  the  phase  points  of 
adjacent  lenslets.  Likewise,  if  the  phase  points  were  all  external  to  the  lenslet  (as  was  shown 
in  Figure  2.1),  there  would  be  an  incorrect  relationship  between  phase  points  and  the  measured 
slope  data  due  to  an  overlap  between  the  neighbor  lenslet  phase  points.  This  arrangement  would 
indicate  that  the  phase  in  one  lenslet  is  influenced  by  the  neighbor,  which  is  not  the  case.  Placing 
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the  phase  points  along  the  edge  in  this  manner  solves  these  two  issues.  The  geometric  impli¬ 
cations  of  how  these  hexagon  lenslets  focus  different  amounts  of  light  per  axis  on  the  sensor 
detector  warrant  additional  inquiry.  The  distance  between  two  neighbor  phase  points  across  a 
row  is  different  than  the  distance  between  two  neighbors  down  a  column.  This  information  must 
be  included  in  the  phase  point  calculation,  or  the  visible  result  will  be  distorted. 

In  the  phase  point  arrangement  depicted  in  Figure  2.4,  a  single  segment  of  the  mirror  can  be 
stored  in  a  rectangular  grid  of  10  rows  and  18  columns.  If  all  six  segments  were  in  a  single 
matrix,  it  can  be  done  in  as  small  of  a  matrix  as  30  rows  and  54  columns. 

2.4.1  Boundary  Conditions 

The  important  consequence  of  the  rectangular  grid  phase  point  arrangement  is  that  the  bound¬ 
ary  correction  algorithm  for  the  circular  aperture  can  now  be  directly  applied  to  the  hexagon. 
Without  this  relationship,  the  formulas  must  be  modified  or  the  slopes  considered  to  lie  along 
the  vectors  that  compose  the  hexagonal  lattice  defined  by  the  center  points  of  the  61  WFS.  The 
actual  formulation  of  the  equations  depends  on  how  the  aperture  boundary  crosses  the  phase 
points.  Every  possible  phase  point  arrangement  should  be  considered;  in  a  rectangular  grid  ar¬ 
rangement  of  four  phase  points,  there  are  16  possible  combinations  of  which  phase  points  are 
inside  and  which  are  outside  of  the  aperture  boundary;  since  each  of  these  is  handled  in  exactly 
the  same  way  as  they  would  for  the  circular  aperture,  no  modification  is  required. 
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CHAPTER  3: 

Fast  Fourier  Transform  Modeling 


The  FFT2  is  an  implementation  of  the  DFT  that  is  computationally  efficient.  In  a  two-dimensional 
grid  of  (j)[m,n],  the  two-dimensional  DFT  provides  the  spatial  frequency  content  of  the  wave- 
front,  showing  that  the  dominant  modes  of  the  wavefront  are  lower  frequencies  [13].  Therefore, 
in  practice  only  the  lower  order  modes  are  needed  to  correct  the  largest  error  terms. 

Breaking  the  wavefront  down  into  spatial  modes  is  well  suited  for  the  actuators.  On  the  SMT, 
the  actuators  are  placed  into  a  triangle  mesh  along  the  structural  supports  on  the  rear  surface  of 
the  mirror.  This  truss  lattice  is  independent  from  the  hexagon  lattice  of  the  WFS.  This  geometric 
relationship  between  the  phase  points  and  the  actuators  varies  across  the  entire  segment  as  seen 
in  Figures  3.1,  C.l  and  C.2.  There  are  two  methods  for  dealing  with  these  different  lattices. 
The  first  is  a  coordinate  transformation.  The  second  method  is  working  in  the  spatial  frequency 
domain. 

In  the  traditional  DM  controller,  the  phase  points  and  actuators  are  co-located.  This  makes  for 
a  simple  controller  because  the  wavefront  phase  error  is  known  at  the  locations  of  the  actuators. 
Thus,  a  simple  actuation  change  causes  a  wavefront  change,  with  the  goal  of  being  an  improve¬ 
ment  in  the  resulting  image.  For  the  SMT  geometry,  this  is  not  the  case.  The  wavefront  phase 
is  sampled  along  the  hexagon  lattice;  the  actuators  are  found  along  the  truss  lattice. 


Figure  3.1:  SMT  wavefront  sensor  and  actuator  topology  for  a  single  segment. 
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y(m,n  -  1  ,t) 


u(m,  n 

y(m,n  +  1  ,t) 

Figure  3.2:  Spring-Mass  system  model  of  a  deformable  mirror. 


y(m  -  1,  n,  t) 


y(m  +  1,  n,  t) 


3.1  Deformable  Mirror  Plant  Model 


The  plant  model  mathematically  represents  the  physical  system,  in  this  case,  the  DM  position. 
To  simplify  the  analysis,  the  mirror  is  thought  of  as  a  system  of  masses  connected  to  one  another 
with  springs.  Each  mass  is  connected  to  four  neighbors,  as  shown  in  Figure  3.2.  In  this  manner, 
the  mirror  surface  forms  a  large  Cartesian  grid  of  mass  points.  Each  point  also  has  a  fifth  spring 
connection  that  is  the  height  displacement  out  of  the  plane  of  the  mirror. 

From  Figure  3.2,  utilizing  the  concept  of  Hooke’s  law,  we  can  relate  the  vertical  displacement 
y(m,n,t),  velocity  y(m,n.t)  and  acceleration  y(m.n.l)  at  location  ( m,n )  and  time  t  to  the  actu¬ 
ators  u(m,n,t)  and  the  neighbor’s  displacement  as 


y(m.  n,t )  =  —  aoy(m,  n,  t )  —  /3oy(m,  n,  t )  +  You(m,n,t ) 


-  (a0,i(y(m,n,t)  -y{m,n  - 

-  (ao-i(y(m,n,t)-y(m,n+l,t))) 

-  (ayo (y{m,n,t) -y(m-l,n,t))) 

-  (a-i,o(,y(m,n,t)  -y(m  +  l,n,t))) 

-  —y{m,n  —  1,?))) 

-  (A),-t (y(m,n,t)  -y(m,n+l,t))) 

~  -y(m-l,n,t))) 

~  03—t.o (y(m,n,t)  -y(m+l,n,t))) 


(3.1) 


where  the  a  coefficients  represent  damping  and  the  /3  coefficients  represent  the  spring  effect 
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from  the  neighbors.  The  negative  signs  are  to  ensure  the  positive  a  and  /3  values  result  in  a 
stable  system. 


The  transfer  function  from  the  input  u(m,n,t )  to  the  output  y(m.  /?,  t)  is  defined  as 


G(m,n,s) 


To 

s2  T  os  T  p  ’ 


(3.2) 


and  the  transfer  function  between  the  center  mass  and  its  neighbors  are  of  the  form 


H(m,n,s) 


—a(m,n)s  —  P(m,n) 
s2  +  as  +  p 


(3.3) 


with  m,n  =  —1,0,1.  We  also  define  a  =  a®  +  Oo.i  +  cto,-i  +  QCi,o  +  a-i,o  and  p  =  /3o  +  /3o,i  + 
A)  ,-i  +Pi.  o  T  /3 —  l.o-  Since  (7  and  p  are  always  positive,  this  system  is  stable.  If  cc  and  j6  do 
not  vary  with  the  index,  the  result  is  a  symmetric  output  where  the  energy  dissipates  into  the 
neighbors  evenly.  These  values  affect  the  pole  locations  of  the  denominator,  and  this  determines 
the  response  time  of  the  system.  In  this  manner,  this  abstract  model  can  be  applied  to  a  DM  if 
these  values  can  be  determined  from  either  structural  analysis  or  testing. 


Equation  (3.1)  can  be  expressed  in  the  time  domain  as 


y(m,n,t) 


-  l  t 

52  [h(h,h,t)y(m^h,n-h,t-T)]+g(t)u(m,n,t 

-/i=-t/2=-t 


d  T 


(3.4) 


which  shows  a  triple  convolution  in  space  and  time.  In  the  time-domain  analysis,  all  neighbor 
actuators  influence  each  other.  The  solution  must  be  obtained  by  concurrently  solving  the  entire 
set  of  equations  for  the  entire  mirror  surface.  Equation  (3.4)  can  be  written  in  the  Laplace 
domain  as 


l  t 

Y(m,n,s)=  52  51  [H(li,h,s)Y(m  —  li,n  —  l2,s)]+G(m,n,s)U(m,n,s)  (3.5) 

h=- t/2=-i 

where  G(m,n,s)  and  H(m.n.s)  are  given  in  Equations  (3.2)  and  (3.3).  For  every  s,  define  the 
two-dimensional  NxN  FFT2  as 
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(3.6) 


r(khk2,s)=  £  £  Y(m,n, 

m= 0 n=0 


and 


H(khk2ls)  H(m,n,s)ezfL^m+k^  (3.7) 

m= 0  /z=0 

with  k\,k2  —  0, . . .  ,N  —  l.  Then  the  discrete  spatial  convolution  in  Equation  (3.5)  becomes 


Y (ki ,k2,s)  —  H(ki,  k2,  s)T  {k\ ,k2,s)  +  G(s)U (k\  ,k2,s).  (3.8) 

This  combination  of  transforms  greatly  simplifies  a  three-dimensional  convolution  into  multi¬ 
plication  in  the  frequency  domain.  In  the  frequency  domain,  each  frequency  is  treated  indepen¬ 
dently  from  the  others;  the  equations  can  be  solved  individually  with  much  less  computation 
required. 

3.2  Deformable  Mirror  Plant  Model  in  State-Space 


Feedback 
Digital  Filter 
F<lf.s) 


Figure  3.3:  Simple  block  diagram  for  feed-forward  and  feedback  transfer  function. 

The  purpose  of  the  mirror  controller  is  to  determine  the  actuator  command  values  at  each  point 
( m,n )  across  the  entire  mirror,  u(m,n,t)  to  remove  as  much  wavefront  error  as  possible.  In  an 
active-optics  system,  the  actuators  might  be  moved  once  during  a  calibration  such  that  the  mirror 
form  is  as  accurate  as  possible,  correcting  the  radius  of  curvature  to  match  designed  radius.  In 
an  adaptive-optics  system,  the  actuators  are  used  as  a  function  of  time  to  correct  for  the  varying 
wavefront  of  fight.  It  is  possible  to  attempt  to  do  both  of  these  corrections  concurrently,  but 
there  would  be  a  loss  of  control  authority  in  that  some  actuators  will  have  a  limited  range  of 
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motion  since  the  actuators  are  already  using  their  stroke  length  for  correcting  the  radius  of  the 
mirror  surface. 

In  order  to  design  the  controller,  first  a  mathematical  model  in  terms  of  the  state-space  equations 
is  derived  for  the  simple  block  diagram  of  transfer  functions  shown  in  Figure  3.3.  The  general 
relationship  between  a  transfer  function  and  its  state-space  equations  is 

G{s)  =  I^f\=C(sI-A)-1B  +  D  (3.9) 

B>o{s) 

where  Nq(s )  is  the  numerator  polynomial  and  Dq(s)  is  the  denominator  polynomial  of  the  trans¬ 
fer  function,  and  A,  B,  C,  D  are  the  matrices  that  relate  the  inputs,  outputs,  and  states.  In  the 
following,  the  D  matrix  is  removed  because  it  is  not  used  for  the  model. 

For  the  feed-forward  and  feedback  blocks,  the  state-space  equations  can  be  derived  as 


Zjk  =  A0zk  +  b0uk, 

(3.10a) 

vk  =  CO  Zk, 

(3.10b) 

Sjc  =  Ask  +  byk, 

(3.10c) 

wk  =  cksk, 

(3.10d) 

yk  =  Vk+wk. 

(3.11) 

where  zk  and  sk  are  the  internal  states  of  the  systems  in  Figure  3.3.  In  Equation  (3.10),  Aq 
relates  the  neighbor  masses  to  one  another,  bo  relates  the  actuators  to  the  masses  and  co  relates 
the  internal  state  to  the  output.  Likewise,  the  coefficients  A,  b  and  ck  do  the  same  for  the 
internal  state  that  varies  in  spatial  frequency.  The  k  subscript  indicates  that  these  equations  and 
coefficients  must  be  applied  to  each  spatial  frequency  vector  <  k\ ,  >  individually  in  actual 

implementation.  These  equations  can  be  combined  to  a  single  state-space  model  to  get 
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(3.12a) 


*k  = 


*k  = 


-Ik  + 


bo 

0 

=8 


Ukl 


(3.12b) 


(3.12c) 


and 


y, t 


(3.13) 


where  fk,  g  and  hk  are  the  equivalent  A,  B  and  C  general  matrices,  respectively.  The  zero  vector, 
0,  indicates  its  row  dimension  will  match  Ao  to  correctly  fill  the  matrix. 

This  state-space  model  is  at  the  basis  of  the  overall  controller  we  will  develop  in  the  next  chapter. 
What  is  significant  about  this  approach  is  that  the  models  represented  by  Equations  (3.12)  and 
(3.13)  are  all  decoupled  and  can  be  individually  controlled. 
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CHAPTER  4: 

Implementation 


In  this  chapter,  we  test  the  estimation  and  control  techniques  developed  using  simulations  in 
MATLAB  and  Simulink  software.  Example  phase  data  sets  were  generated  by  algorithms  de¬ 
veloped  by  Melissa  Corley  [14]  while  with  the  Spacecraft  Research  and  Design  Center  (SRDC) 
at  NPS.  Her  code  has  excerpts  from  Chris  Wilcox  on  Spatial  Light  Modulators  (SLM)  [15]. 
The  simulated  phase  data  size  of  the  aperture  was  1024x1024  phase  points,  much  larger  than 
existing  adaptive-optics  systems  at  the  SRDC.  This  large  dimension  was  chosen  to  show  the 
scalability  of  the  algorithm. 

4.1  Wavefront  Reconstruction 

The  circular  aperture  solution  forms  the  basis  for  solving  a  number  of  different  geometries. 
Poyneer  outlined  the  general  theory  of  how  to  accomplish  this  [6].  In  this  section,  a  particular 
algorithm  for  this  technique  is  developed. 

One  of  the  difficulties  of  applying  a  Fourier-based  reconstruction  algorithm  is  the  fact  that 
apertures  are  generally  not  square  but  circular  or  polygonal,  which  does  not  fit  with  the  square 
grid  of  a  FFT2.  In  order  to  overcome  this  difficulty,  we  need  to  make  some  assumptions  on  the 
data  outside  the  aperture  and  within  the  grid  of  the  FFT2.  In  general,  we  assume  the  extra  values 
to  be  zero.  Particular  care  has  to  be  taken  at  the  boundary  of  the  aperture  since  the  approach 
developed  by  Poyneer  leads  to  algorithms  which  are  fairly  sensitive  to  boundary  conditions  [6] . 

Apertures  with  various  characteristics  of  vignetting,  such  as  mechanical  obstructions  to  the  op¬ 
tical  pathway,  can  also  be  included  in  the  final  solution.  This  can  be  useful  for  space-based 
telescope  apertures  since  the  secondary  mirror  usually  has  mechanical  vignetting  from  its  struc¬ 
tural  support  that  holds  the  mirror  in  place. 

To  implement  this  algorithm,  the  entire  grid  of  data  is  iterated  to  identify  each  group  of  four 
neighbor  phase  points.  The  only  groups  of  consequence  are  along  the  boundary.  The  algorithm 
builds  a  list  of  boundary  slopes  u  that  needs  to  be  solved.  If  another  new  unknown  is  found 
along  the  boundary,  it  is  appended  sequentially  in  a  list.  If  it  is  already  in  the  list,  the  algorithm 
uses  the  number  previously  assigned.  For  the  square,  any  slopes  that  are  completely  internal  to 
the  aperture  need  to  be  used  in  computing  the  constant  for  c.  After  the  entire  grid  is  processed, 
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Figure  4.1:  Square  of  phase  points  for  generating  equations. 

all  of  the  equations  are  identified.  The  matrix  M  in  Equation  (2.9)  is  quickly  produced.  All  of 
the  required  information  for  Equation  (2.9)  is  available  to  solve  the  boundary  problem.  This 
design  allows  for  a  variety  of  apertures  to  be  solved,  including  the  hexagon. 

In  Figure  4.1,  the  black  dashed  line  shows  the  aperture  boundary.  The  blue  phase  point  dot  is 
outside  of  the  aperture,  while  the  black  dots  are  inside.  The  red  dotted  line  indicates  the  mesh 
equation  direction.  The  unknowns  and  internal  slopes  are  identified.  In  this  example,  the  entry 
in  c  in  Equation  (2.9)  is  equal  to  sy[j,k+  1]  —  sx[j  +  1  ,k\. 

After  correctly  setting  up  Equation  (2.9)  variables  M  and  c,  u  is  easily  solved.  The  values  u 
contains  apply  to  specific  sx[j.k]  or  sy[j.  k]  that  are  identified  by  the  unknowns.  Each  u  value  is 
iterated  and  then  applied  to  the  slopes  in  Equation  (2.8).  The  results  are  matricies  of  ix[j,  k]  and 
iy[j,k]  which  are  used  in  the  reconstruction  algorithm.  This  technique  is  applied  to  the  large 
example  data  set  shown  in  Figure  1.2.  The  reconstructed  wavefront  is  shown  in  Figure  4.2.  The 
reconstruction  looks  comparable  to  the  original  phase  data,  but  there  are  some  artifacts  in  the 
external  aperture  area. 


Sy\j,k+1\ 


Figure  4.2:  (a)  Original  wavefront;  (b)  Reconstructed  wavefront. 
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The  difference  in  the  reconstructed  wavefront  from  the  original  (known)  wavefront  is  computed 
and  shown  in  Figure  4.3.  This  shows  that  the  larger  errors  are  at  the  boundary  as  expected 
or  at  the  external  edge  of  the  boundary.  This  is  a  consequence  of  the  filtering,  discarding  the 
imaginary  components  and  round-off  errors  in  floating  point  operations.  Since  the  algorithm  is 
insensitive  to  the  mean  value,  the  comparison  is  made  between  the  reconstructed  data  and  the 
mean- subtracted  original  data. 


Figure  4.3:  Error  in  the  wavefront  reconstruction. 


To  better  compare  the  results,  the  reconstructed  data  is  further  processed  by  zeroing  out  any 
entries  outside  of  the  known  aperture.  In  the  actual  implementation,  these  values  are  not  used 
or  needed.  As  can  be  seen  in  Figure  4.4,  the  estimate  data  closely  matches  the  original  data. 
These  additional  steps  are  not  necessary  in  actual  implementation  and  are  only  included  here  to 
help  the  reader  see  that  the  reconstruction  was  successful. 

4.2  Fourier  Transform  Modeling 

The  DM  is  modelled  using  the  discrete  spatial  Fourier  transform  (DSFT)  to  decompose  the 
spatial  modes  of  the  mirror.  Mathematically,  this  is  done  using  rectangular  matrices.  Matrix 
sizes  are  padded  with  zeros  to  a  length  that  is  a  multiple  of  two.  This  is  done  so  that  the  FFT2 
algorithm  can  execute  efficiently.  Previous  work  has  been  done  using  64x64  matrices  [6].  In 
this  thesis,  the  matrices  used  were  256x256  in  size  to  show  that  the  algorithm  is  scalable  to  the 
larger  apertures  that  are  expected  in  the  future. 


23 


Figure  4.4:  Cleaned  up  estimated  wavefront  reconstruction. 
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Figure  4.5:  Plant  model  for  deformable  mirror. 

The  plant  implementation  can  be  seen  in  Figure  4.5.  In  order  to  make  the  implementation  more 
efficient,  the  analog  transfer  functions  in  s  are  implemented  in  discrete  time  in  z  by  converting 
from  continuous  time  to  discrete  time  using  a  zero-order  hold  (ZOH)  over  an  assigned  sampling 
period  Ts.  The  sampling  time  can  be  adjusted  to  match  the  physical  system  response.  In  order 
to  avoid  algebraic  loops,  we  implement  the  transfer  functions  to  be  strictly  proper  by  adding  a 
time  delay  1 . 

The  forward  gain  in  Equation  (3.2)  is  implemented  as  a  scalar  multiplier  of  each  spatial  fre¬ 
quency.  To  implement  Equation  (3.3),  two  matrices  are  created,  ha  and  hp.  Each  matrix  size  is 
set  to  be  the  same  as  the  size  of  the  aperture  and  has  the  form 


24 


0  -«o,i  0  ...  0 

-aio  0  0  ...  0 

ha(m,n,s) 

o  o  o  0 

-a_po  0  0  ...  0 

After  creating  these  matrices,  their  two-dimensional  FFTs  are  computed.  The  two  matrices  ha 
and  h  a  are  needed  for  the  Laplace  s  term  and  the  constant  term.  These  matrices  are  constructed 
with  periodicity.  The  term  that  relates  the  mirror  mass  to  its  northern  neighbor  is  wrapped 
around  to  the  bottom  row  of  the  first  column;  likewise,  the  term  that  relates  the  mirror  mass 
to  its  western  neighbor  is  wrapped  around  to  the  rightmost  entry  in  the  first  row.  The  DSFT 
is  applied,  and  the  results  are  identified  as  matrices  Ha  and  Hp .  The  three  matrices  are  used 
in  conjunction  with  the  digital  filters  as  an  element-wise  matrix  multiplier,  which  implements 
Equation  (3.8)  efficiently  as  compared  to  Equation  (3.4).  These  matrices  are  not  converted  to 
discrete  time,  as  the  discrete  filters  implement  the  continuous-time  equation. 

The  purpose  of  the  controller  is  to  meet  the  design  goal  of  Equation  (4.2)  given  that  e(m,n,t)  is 
the  error  of  the  mirror  state  minus  the  sensed  wavefront.  The  ability  of  the  system  to  keep  the 
error  below  this  threshold  depends  upon  how  rapidly  the  mirror  can  respond  to  the  wavefront’s 
constantly  evolving  phase  progression  in  time.  This  requirement  is  expressed  by 


e(m,n,t)  =  y(m,n,t)  —  (j)(m,n,t)  (4.2a) 

e{m,n,t) \  <  £,Vt  >  to.  (4.2b) 


The  last  step  is  to  add  an  integrator  state  to  the  final  output  of  the  controller.  An  integrator 
controller  has  additional  robustness  compared  to  classical  controllers.  The  integral  controller 
is  well  suited  for  tracking  of  constant  reference  inputs,  which  in  the  case  of  slowly- varying 
wavefront  mode  coefficients  is  an  adequate  approximation.  The  resulting  state-space  equations 


are 


0  0  Uk 

§  fk  xk 


=fk  =s 


(4.3) 


25 


and 


yk  = 


=hk 


Uk 

■Ik 


(4.4) 


For  each  spatial  frequency,  the  matrix  dimensions  for  fk,  g  and  hk  are  5x5,  5x1  and  1x5,  re¬ 
spectively.  However,  this  fifth  order  system  is  not  of  full  rank.  An  uncontrollable  mode  exists 
because  of  a  pole/zero  cancellation  in  the  transfer  function.  Referring  to  Figure  3.3,  we  obtain 
the  transfer  function  which  shows  this  cancellation: 

Yk(s)  No(s)  1  No(s)  D(s)  No(s) 

Uk(s)  D(s)i_ D(s)  D(s)-Bk(s)  D(s)-Bk(s)' 

D(s) 

The  pole-zero  cancellation  is  stable  and  does  not  affect  the  stability  of  the  system.  The  state- 
space  equivalent  matrices  can  be  computed  from  the  transfer  functions  which  directly  define 
the  digital  filters.  The  matrices  fk,  g  and  hk  can  all  be  defined  from  Equation  (4.5)  using  Equa¬ 
tion  (3.9).  Then  matrices  fk ,  g  and  hk  can  be  defined  from  Equations  (4.3)  and  (4.4).  At  this 
point,  the  controller  is  well  defined,  but  an  observer  needs  to  be  added.  The  observer  estimates 
the  mirror  height  displacement  position  at  each  mass  position  since  in  the  real  system  the  true 
mirror  height  displacement  is  not  known.  The  control  signal  uk  is  computed  by  combining  a 
state  estimator  and  state  feedback  as 


xk  =  fkxk  +  guk  +  Kk(yk  -  hkxk)  (4.6) 

and 

uk  =  ~Lk%  (4.7) 

with  Kk  such  that  det (si  —  fk  +  Kkhk )  =  r(.s')  and  Lk  such  that  det {si  —  fk  +  gLk )  =  A(s).  Both 
of  the  polynomials  F(.v)  and  A(.v)  should  be  exactly  the  same  for  all  k  in  <  k \  Mi  >■  In  order 
to  accomplish  this,  the  values  for  Kk  and  Lk  must  be  solved  for  each  k  by  selecting  the  pole 
placement  locations.  These  placements  can  be  made  to  simulate  the  expected  real-world  system 
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N<k  .2) 

Digital  Filters 


Figure  4.6:  Controller  model  for  deformable  mirror, 
with  additional  information  from  the  SMT. 

All  of  the  system  developments  up  to  now  lead  to  the  creation  of  the  system  shown  in  Figure  4.6. 
The  information  needed  to  produce  the  combination  feed-forward  and  feedback  system  that  can 
control  the  DM  is  highlighted  by  this  model.  The  observer  is  made  to  follow  the  rule  of  thumb 
that  it  converge  to  the  actual  state  at  four  to  ten  times  faster  than  the  state  changes  in  the  pole 
selection.  The  number  of  poles  placed  in  both  the  controller  and  observer  equals  the  order 
of  the  system,  which  in  this  case  is  three.  In  the  figure,  the  W  and  V  matrices  represent  the 
polynomial  coefficients  of  the  transfer  function  numerator  terms.  There  are  three  each,  while 
the  denominator  is  a  fourth  order  polynomial  that  is  the  same  for  all  k. 

Combining  Figures  4.5  and  4.6,  we  get  the  resulting  complete  simulation  shown  in  Figure  4.7. 
The  “phi”  data  input  into  the  system  shown  in  Figure  4.7  is  from  a  phase  generation  program  run 
externally  before  the  simulation  [14].  Ideally,  in  the  real  system,  the  wavefront  reconstruction 
algorithm  discussed  in  this  thesis  can  be  used  to  generate  the  estimated  phase. 

The  simulation  is  configured  with  demonstration  parameters  that  exhibit  the  dynamics  expected 
from  a  DM.  The  resulting  system  is  stable,  and  the  time- varying  wavefront  is  fed  into  the  system 
controller  which  provides  the  actuator  displacements  to  the  mirror. 

In  Figure  4.8,  the  first  25  wavefront  updates  are  shown,  and  the  single  pixel  maximum  error  is 
plotted.  The  controller  is  adjusting  the  actuators  to  correct  for  the  error  as  shown  in  the  plot.  The 
wavefront  is  sampled  at  a  lower  rate  than  the  actuators  can  be  moved,  so  the  saw-tooth  pattern 
shows  that  the  error  is  approaching  zero.  When  the  new  wavefront  is  provided,  the  maximum 
pixel  error  has  a  discontinuous  jump,  and  the  pattern  keeps  repeating.  The  percentage  is  cal¬ 
culated  based  on  the  absolute  maximum  phase  value  over  the  entire  aperture.  This  information 
can  be  treated  as  a  figure  of  merit  to  show  that  the  controller  is  working  across  the  large  set  of 
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Figure  4.7:  Simulation  model  for  deformable  mirror  controller. 


Maximum  Pixel  Error  Percentage 


Figure  4.8:  Maximum  pixel  error  across  the  entire  aperture  for  25  wavefront  updates. 


data. 

In  Figure  4.9,  the  pixel  error  is  broken  up  into  the  difference  between  the  maximum  and  the 
minimum.  This  shows  the  error  envelope.  The  negative  error  is  much  larger  initially  from  the 
large  phase  changes  in  the  blue  regions  shown  in  Figure  1.2.  The  error  on  both  sides  now  has 
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Bounded  Pixel  Error  Percentage 


Figure  4.9:  Maximum  error  extents  on  either  side  of  the  wavefront  across  the  entire  aperture 
for  25  wavefront  updates. 

the  saw-tooth  pattern.  At  the  beginning  of  the  plot,  the  lower  side  error  makes  large  jumps  from 
larger  changes  in  the  example  wavefront.  Despite  the  large  changes,  the  system  responds  by 
a  much  sharper  rate  of  error  decrease.  This  is  an  interesting  point.  The  graph  results  are  not 
tracking  the  same  individual  pixels  throughout  the  simulation  but  rather  the  specific  pixels  with 
maximum  error  extents. 

In  the  spatial  frequency  domain,  the  majority  of  the  wavefront  energy  is  concentrated  in  the 
lower  frequencies.  This  detail  can  be  exploited  further  to  minimize  the  required  computing  of 
all  spatial  frequencies  to  just  the  range  of  interest.  If  the  system  can  tolerate  the  higher  fre¬ 
quency  content  being  uncontrolled,  the  computing  requirements  will  decrease.  In  this  manner, 
extremely  large  apertures  can  be  processed  without  higher  end  computing  resources.  This  can 
also  be  exploited  on  a  space-based  telescope,  where  the  processors  lag  ground-based  perfor¬ 
mance  by  nearly  a  decade. 
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CHAPTER  5: 
Conclusions 


5.1  Summary 

The  main  focus  of  this  thesis  is  the  application  of  the  SMT  to  adaptive  optics.  In  particular 
it  addressed  two  important  issues:  wavefront  reconstruction  and  feedback  control.  Wavefront 
reconstruction  is  developed  based  on  the  fundamental  geometry  of  the  data  and  the  grid  points 
where  the  wavefront  phase  is  measured.  Interpolating  between  these  points,  we  compute  the 
estimate  of  the  wavefront.  Techniques  vary  by  computational  complexity.  While  ground-based 
telescopes  can  have  large  computing  power,  space-based  telescopes  usually  have  other  limit¬ 
ing  factors  of  size,  weight,  power  and  electronics  reliability  to  consider.  While  the  industrial 
community  is  making  larger  apertures  with  increased  actuator  density,  the  algorithm’s  compu¬ 
tational  requirements  should  not  scale  faster  than  the  computing  capabilities  of  the  platform  in 
future  systems. 

The  algorithms  presented  in  this  research  are  based  on  the  DFT,  which  is  computationally  ef¬ 
ficient  and  offers  results  comparable  to  other  techniques  such  as  the  Zernike  Polynomials.  The 
wavefront  can  be  reconstructed  for  non-noisy  data  or  with  a  least-squares  fit  for  noisy  data. 

Similarly,  the  control  is  based  on  a  simple  model  of  the  DM  as  a  mesh  of  masses  connected  by 
springs.  The  mirror  has  actuators  along  the  normal  of  the  rear  surface  that  adjust  the  heights. 
The  model  is  developed  from  basic  theory  and  has  not  yet  been  adopted  yet  to  the  SMT.  By 
assuming  linearity  and  stationarity  in  the  time  and  spatial  domains,  we  can  take  full  advantage 
of  a  decomposition  in  terms  of  spatial  frequencies  using  the  FFT2.  In  this  way,  the  highly 
connected  time-space  domain  model  is  reduced  to  a  set  of  independent  models  in  the  frequency 
domain  which  can  be  easily  controlled. 


5.2  Additional  Work 

The  wavefront  reconstruction  was  accomplished  for  the  single  hexagon  panel.  Additional  work 
could  apply  this  to  the  six  panels  of  the  SMT.  The  interesting  contribution  of  this  would  be  to 
determine  if  the  panels  should  be  solved  independently  or  as  a  larger  more  complicated  aperture 
to  see  under  which  circumstances  would  the  error  would  improve. 
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In  addition,  there  are  a  number  of  tests  to  run  on  the  actual  telescope  hardware.  The  effect  of 
mechanical  vignetting  of  the  SMT  in  wavefront  reconstruction  is  not  known.  The  relationship 
of  the  outer  edge  of  the  hexagon  aperture  to  the  mirror’s  flat  edges  to  the  lenslet  jagged  hexagon 
edges,  seen  in  Figure  3.1,  is  not  known.  The  lenslets  may  all  be  contained  inside  of  the  mirror 
edge  or  extend  pass  the  edge.  If  there  is  overlap,  it  should  be  studied  to  determine  the  effect  on 
wavefront  reconstruction. 

There  has  been  previous  work  done  by  Puschel  and  Rotteler  on  the  discrete  Triangle  Trans¬ 
form  (DTT)  [16].  This  work  could  be  applied  to  the  truss  lattice  of  the  actuators.  Additionally, 
Vince  and  Zheng  have  worked  on  taking  the  DFT  of  a  hexagon  lattice  [17].  This  work  could  be 
applied  to  the  wavefront  sensor  lattice.  Future  work  could  compare  which  method  allows  for 
an  efficient  transformation  between  lattices  and  solve  the  varying  relationships  for  the  actuators 
to  wavefront  sensors  in  the  lattices.  There  could  be  interesting  consequences  of  the  DTT  ap¬ 
proach  given  that  it  is  similar  to  the  discrete  Cosine  Transform  (DCT-III)  and  is  less  sensitive  to 
boundary  conditions.  A  study  of  how  the  data  processing  requirements  vary  from  this  technique 
would  be  beneficial. 

An  interesting  consequence  of  the  hexagonal  lenslets  is  that  they  are  not  symmetric  about  the 
two-slope  axis.  This  could  have  interesting  effects  on  the  wavefront  reconstruction  and  sec¬ 
ondary  camera,  a  topic  worth  exploring  with  additional  research. 
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APPENDIX  A: 

Equation  Set  for  Circular  Aperture 


The  following  equations  represent  a  solution  to  the  boundary  conditions  represented  by  Equa¬ 
tions  (2.8)  for  Figure  2.3: 


U2  —  “1  =  0 
U4+Sy[ 5, 1]  -  sx[ 4, 2]  —  u3  =  0 
“7  -sx[6,2]  -m6  =  0 
“8  +  u9  -  £*[ 7, 3]  -  [7, 2]  =  0 
M  | o  —  “ I  I  -5y[8,3]  =  0 
yv[8,5]  +  ni3  -Mi4-5V[8,5]  =  0 

“15  —  “16  =  0 
“17  —  “18  =  0 

yv[5, 8]  +  w19  -  u2o  -  5y[5, 8]  =  0 
■Sjc[3,  8]  +  U22  —  “23  =  0 
Sx[2,7]+Sy[3,7]-U24~U25  =0 

“27 +  ^[2, 6]  “26  =  0 

U30  +  Sy[2,4\  - sx[l,5]  -  U29  =  o 

— «32  —  «31  =  0 


“3  Ty[3,  2]  m2  =0 

“5+ “6 -^[5,2]  -sy[5,l]  =0 

—  —  U7=0 

—  MiO  —  Ug  =  0 
“ii  +“12 -s*[8,5]  -5v[8,4]  =0 
«14  —  U 15  —  5y[8,6]  =  0 
sx[l,l]  +  ui6-un-sy[l,l]  =0 
5x[6,8]  +M18  —  “19  =0 
sx[ 4, 8]  +  .s-y [5,8]  -  m2i  -  i<22  =  0 

“24  +  “23  =  0 
“26  +  “25  =  0 

yY[l,5]+5v[2,5]  -m27-“28  =0 
“31  +  Sv[2, 3]  —  w30  =  0 
Ml  +5-y[3,2]  —  sx[2,3]  —  m32  =0 


These  are  shown  for  completeness  to  illustrate  the  example.  These  equations  are  created  by 
writing  the  mesh  equation  for  each  square.  The  path  followed  along  each  square  is  top,  right, 
bottom,  left.  If  the  slope  arrow  points  in  the  opposite  direction  of  this  path,  it  is  distinguished 
by  negative  sign. 
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APPENDIX  B: 
MATLAB  code 


B.l  ComputeSlopes.m 

%  This  file  implements  an  algorithm  for  solving  the  boundary  slopes 
%  and  then  displays  the  results  for  a  sample  set  of  data. 

%  Useful  output  variables: 

%  unknowns  —  structure  of  the  unknown  slopes,  fields  are:  axis,  j,  k 
%  equations  —  array  of  the  unknowns  for  each  equation  (2  per  eqn) 

%  fields  are:  axis,  j,  k,  sign 

%  knowns  —  array  of  known  slopes  needed  for  each  equation  (2  per  row) 
%  This  is  useful  for  dynamic,  time  varying  data 

%  Just  iterate  over  this,  sum  the  values  per  row  for  new  c  vector 
%  fields  are:  axis,  j,  k,  sign 

%  M,  Minv  —  matrices  that  correspond  to  equations. 

%  code  by  Travis  Axtell 

close  all; 
clear  all; 
clc; 

%  load  phi  variable 
load  data; 

%  pad  zeros  on  all  sides  to  guarantee  an  edge 
[sizel,  size2 ]  =  size(phi); 
phi2  =  zeros ( sizel+2 , size2+2 ) ; 
phi2 (2 :  (sizel  +  1) , 2 :  (size2  +  l) )  =  phi; 
phi  =  phi2; 
clear  phi2; 

%  Find  the  aperture  of  data.  This  code  depends  on  the  actual  data. 

%  For  this  example  data,  the  edge  is  entirely  surrounded  by  0. 

%  This  should  be  confirmed  before  applying  new  data  to  algorithm, 
aperture  =  phi  ~=  0; 
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53 

54 

55 
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57 
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%  Find  the  slopes  from  the  example  phase  data. 

%  Real  world  data  would  just  be  provided  without  knowing  the  phase, 
sx  =  [phi ( : ,  2:sizel+2),  phi ( : ,  1)]  —  phi; 
sy  =  [phi  (2 : sizel+2 ,  :);  phi(l,  :)]  —  phi; 

%  Empty  variables  to  begin  the  code. 

unknowns  =  [];  %  list  of  structures  that  identify  the  unknown  slopes 
knowns  =  [];  %  list  of  structurs  that  identify  the  known  slopes 
sums  =  [];  %  vector  of  sums,  c,  to  solve  for  u. 
equations  =  [];  %  list  of  equations  using  the  unknowns 

%  Paceholder  for  zeros  as  knowns . 

kzero  =  struct ( 1  axis  * ,  'zero’,  'j',  j,  'k',  k,  'sign',  1) ; 

%  iterate  over  each  square  of  values  from  the  aperture 
for  j  =  l:sizel+l, 

for  k  =  l:sizel+l, 

%  binary  value  for  the  arrangement  of  the  4  values  are  as  shown: 

%  1  2 
%  4  8 

value  =  aperture ( j , k)  +  2*aperture  (  j , k+1 )  +  ... 

4*aperture  (  j  +  1 , k)  +  8*aperture  (  j  +  1 , k+1 ) ; 

%  This  switch  statement  is  not  truly  required  but  helps  to  organize 
%  the  code . 
switch  value 

case  {0,  15  } 

%  the  block  is  all  internal  or  external  and  no 

%  processing  is  required. 

continue; 

case  {  1,  2,  4,  8  } 

%  single  point  on  aperture. 

%  Each  if  statement  below  identifies  the  two  unknowns  and 
%  their  associated  mesh  equation  sign, 
sum  =  0; 

%  kl,  k2  are  placeholders  since  they  are  zeros, 
kl  =  kzero; 
k2  =  kzero; 
if  value  ==  1, 

ul  =  struct (* axis 1 ,  'x',  'j',  j,  'k',  k,  'sign1,  1) ; 
u2  =  struct  (' axis ' ,  ’y',  1 j ' ,  j,  *k',  k,  'sign',  —1); 


36 
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91 

92 
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94 
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101 

102 
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110 
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case 


elseif 

value  ==  2, 

ul 

=  struct ( ' axis ' , 

'x'. 

j,  - k 1 

f 

r 

k,  ' 

sign ' , 

1 

)  ; 

u2 

=  struct ( ' axis  ' , 

'y\ 

'  j 

j,  ' k 1 

i 

r 

k+1 , 

' sign ' 

r 

i) ; 

elseif 

value  ==  4, 

ul 

=  struct ( ' axis ' , 

'y'. 

'  j'. 

j,  •  k 1 

i 

r 

k,  ' 

sign ' , 

- 

i) ; 

u2 

=  struct ( ' axis ' , 

'X', 

'  j'. 

j+i,  ' 

'  k ' 

1  ,  k. 

' sign ' 

r 

-i)  ; 

elseif 

value  ==  8, 

ul 

=  struct ( ' axis ' , 

'X', 

'  j'. 

j+i,  ' 

'  k  ’ 

’  ,  k. 

' sign ' 

r 

-i)  ; 

u2 

=  struct ( ' axis ' , 

'Y\ 

'  j'. 

j,  ' k 1 

f 

r 

k+1. 

' sign ' 

r 

i) ; 

end 


%  checkadd  confirms  the  unknown  is  new  before  adding  it  to 
%  the  unknowns  list,  or  if  it  is  not  new,  redefines  the  ul 
%  to  have  the  correct  unknown  number  for  generating  the 
%  equations. 

[unknowns,  ul]  =  checkadd (unknowns ,  ul)  ; 

[unknowns,  u2]  =  checkadd (unknowns ,  u2)  ; 


{  3,  5,  10,  12  } 

%  two  points  on  aperture, 
sum  =  0; 

%  Each  if  statement  below  identifies  the  two  unknowns  and 
%  their  associated  mesh  equation  sign, 
if  value  ==  3, 


ul 

=  struct ( ' axis  1 

r 

'y'r 

'  j\ 

j,  ' k ■ 

1 

r 

k,  'sign',  -1); 

u2 

=  struct ( ' axis  1 

r 

'y'. 

'  j', 

j,  •  k ■ 

1 

t 

k+1 ,  1  sign ' ,  1 ) ; 

sum 

l  =  sx ( j,  k)  ; 

kl 

=  struct ( ' axis ' 

r 

'X', 

'  j', 

j,  ' k ■ 

1 

r 

k,  ' sign ' ,  1) ; 

k2 

=  kzero; 

elseif 

value  ==  5, 

ul 

=  struct  (  ' axis  1 

r 

'X', 

'  j', 

j,  ’ k ■ 

• 

r 

k,  ' sign ' ,  1) ; 

u2 

=  struct  (  ' axis  1 

r 

'X', 

'  j\ 

j+i,  ' 

'  k ' 

1 ,  k,  ' sign ' ,  — 1) 

sum 

i  =  — sy ( j , k) ;  % 

negative  for  mesh 

equation 

kl 

=  struct  (  ' axis  1 

r 

'y'r 

'  j', 

j,  -k' 

1 

r 

k,  'sign',  -1); 

k2 

=  kzero; 

elseif 

value  ==  10, 

ul 

=  struct ( ' axis  1 

r 

'x'. 

'  j', 

j,  ' k ■ 

1 

t 

k,  ' sign ' ,  1) ; 

u2 

=  struct  (  ' axis ' 

r 

'X', 

'  j', 

j+i,  ' 

1  k ' 

!,  k,  'sign',  -1) 

sum 

i  =  sy ( j, k+1) ; 

kl 

=  struct  (  ' axis ' 

\ 

f 

’y\ 

'  j', 

j,  -k' 

1 

t 

k+1 ,  ' sign ' ,  1 ) ; 

k2 

=  kzero; 

elseif 

value  ==  12, 
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137 

138 
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141 
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150 
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156 


case 


ul  = 

struct ( ' axis ' , 

'y ', 

'  j'. 

j,  ' 

k'  , 

k,  ' sign ' , 

-i) ; 

u2  = 

struct ( ' axis ' , 

'Y'r 

'  j'. 

j,  ' 

k'  , 

k+1 ,  ' sign ' 

,  D  ; 

sum  : 

=  — sx ( j+1,  k) ;  % 

negative 

for  mesh 

equation 

ll  II 

i—l  CM 

m  2*: 

struct ( ' axis ' , 

kzero ; 

'x'. 

'  j'. 

j  +  1/ 

'k' 

,  k,  'sign' 

,  -i); 

end 

%  checkadd  confirms  the  unknown  is  new  before  adding  it  to 
%  the  unknowns  list,  or  if  it  is  not  new,  redefines  the  ul 
%  to  have  the  correct  unknown  number  for  generating  the 
%  equations. 

[unknowns,  ul]  =  checkadd (unknowns ,  ul)  ; 

[unknowns,  u2]  =  checkadd (unknowns ,  u2) ; 

{  7,  11,  13,  14  } 

%  three  points  on  aperture, 
sum  =  0; 

%  Each  if  statement  below  identifies  the  two  unknowns  and 


%  their 

associated  mesh 

equation 

sign . 

if  value 

==  7, 

ul  = 

struct ( ' axis 

r 

'X', 

'  j’ 

,  j  +  1/ 

'k' 

,  k,  'sign',  -1) ; 

u2  = 

struct ( ' axis 

r 

'Y'r 

'  j' 

,  j/  'k 

! 

r 

k+l ,  ' sign ' ,  1 ) ; 

sum 

=  sx ( j , k)  -  sy  (  j 

,k)  ; 

%  signs  for  mesh  equation 

kl  = 

struct ( ' axis 

r 

'X', 

'  j’ 

,  j/  'k 

! 

r 

k,  ' sign ' ,  1) ; 

k2  = 

struct  (  ' axis 

r 

'y'r 

'  j' 

,  j/  'k 

! 

r 

k,  'sign',  -1) ; 

elseif  value  ==  11, 

ul  = 

struct  (  ' axis 

r 

'y'. 

'  j' 

,  j/  'k 

i 

r 

k,  'sign',  -1) ; 

u2  = 

struct  (  ' axis 

r 

'X', 

'  j' 

/  j  +  1. 

'k' 

,  k,  'sign',  -1) ; 

sum 

=  sx ( j,  k)  +  sy ( j 

,  k+l ) 

} 

signs  . 

for 

'  mesh  equation 

kl  = 

struct ( ' axis 

r 

'X', 

'  j' 

,  j,  'k 

I 

r 

k,  ' sign ' ,  1) ; 

k2  = 

struct ( ' axis 

r 

'y'z 

'  j' 

,  j,  'k 

! 

r 

k+l ,  ' sign ' ,  1 ) ; 

elseif  value  ==  13, 

ul  = 

struct ( ' axis 

r 

'X', 

'  j' 

,  j,  'k 

! 

r 

k,  ' sign ' ,  1) ; 

u2  = 

struct ( ' axis 

r 

'y'r 

'  j' 

,  j,  'k 

! 

r 

k+l ,  ' sign ’ ,  1 ) ; 

sum 

=  — sy ( j , k)  - 

sx ( j+1, k) ; 

%  signs 

for  mesh  equation 

kl  = 

struct ( ' axis 

r 

'y'z 

'  j' 

,  j,  'k 

! 

r 

k,  'sign',  -1) ; 

k2  = 

struct ( ' axis 

r 

'X', 

'  j' 

,  j  +  1. 

'k' 

,  k,  'sign1',  -1)  ; 

elseif  value  ==  14, 

ul  = 

struct ( ' axis 

r 

'X', 

'  j' 

,  j,  'k 

1 

r 

k,  ’ sign ’ ,  1) ; 

u2  = 

struct ( ' axis 

? 

t 

'y\ 

'  j' 

,  j,  'k 

i 

r 

k,  'sign',  -1); 

sum 

=  sx (j+1,  k) 

+ 

sy  ( j. 

k+l ) ;  %  signs 

for  mesh  equation 

kl  - 

struct ( ' axis 

r 

'X', 

'  j' 

,  j  +  1. 

'k' 

,  k,  'sign',  -1); 
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k2  -  struct ( 1  axis ' ,  'y',  'j',  j ,  'k',  k+1,  'sign1,  1); 

end 

%  checkadd  confirms  the  unknown  is  new  before  adding  it  to 
%  the  unknowns  list,  or  if  it  is  not  new,  redefines  the  ul 
%  to  have  the  correct  unknown  number  for  generating  the 
%  equations. 

[unknowns,  ul]  =  checkadd (unknowns ,  ul) ; 

[unknowns,  u2]  =  checkadd (unknowns ,  u2)  ; 

otherwise 

%  This  code  should  never  get  executed  since  all  cases  are 
%  handled,  but  is  included  to  verify  that  the  code  executes 
%  correctly. 

%  W:  for  warning 

fprintf('W:  (j,  k)  =  (%d,  %d) ;  value  =  %d\n',  j,  k,  value); 

end 

%  The  lists  are  updated  to  include  the  new  entries 
equations  =  [equations;  ul,  u2]  ; 
sums  =  [  sums;  sum  ] ; 
knowns  =  [  knowns;  kl,  k2  ]; 

end 

end 

%  Develop  the  M  matrix  from  the  equations 
eqnlen  =  length (equations ) ; 

M  =  zeros (eqnlen, eqnlen) ; 

for  m  =  1: eqnlen, 

%  Each  row  of  M  contains  1  equation,  which  includes  2  entries. 

M (m,  equations (m, 1 ). number )  =  equat ions (m, 1 ). sign; 

M (m,  equat ions (m, 2 ). number )  =  equations (m, 2 ). sign; 

end 

%  If  the  Minv  is  needed  to  be  computed,  uncomment  this 
%Minv  =  pinv (M) ; 

%  In  this  case,  the  Minv  was  precomputed.  This  is  how  it  would  be  done  in 
%  actual  implementation  for  faster  results, 
load  Minv. mat; 

%  Solving  for  u  is  now  a  simple  statement. 
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%  For  dyanmic  data,  rebuild  sums  using  knowns  before  this  next  line, 
u  =  Minv*sums; 

%  Now  have  a  solution.  Just  need  to  apply  it  to  the  original  data. 

%  sx  =  ix  +  bx  =>  ix  =  sx  —  bx 

%  ix  internal  slope,  sx  computed  slope,  bx  boundary  slope 

%  ix  is  what  is  needed,  sx  is  what  was  measured,  bx  is  the  error  to  remove 
ix  =  sx; 
iy  =  sy; 


unklen  =  length (unknowns ) ; 
for  m  =  l:unklen, 

%  This  if  statement  determines  which  set  of  data  to  apply  the  solved 
%  unknown  to. 

if  unknowns (m) .axis  ==  '  x', 

ix (unknowns (m) . j ,  unknowns (m) . k)  =  ... 

sx (unknowns (m) . j ,  unknowns (m) . k)  —  u (m) ; 
elseif  unknowns (m) . axis  ==  '  y', 

iy (unknowns (m) . j ,  unknowns (m) . k)  =  ... 

sy (unknowns (m) . j ,  unknowns (m) . k)  —  u (m) ; 

end 

end 

%  The  mean  is  removed  from  the  data  before  its  FFT2  is  taken. 
sXfft  =  f ft2 ( ix— mean (mean ( ix) ))  ; 
sYfft  =  f ft2 (iy— mean (mean (iy))); 


N  =  sizel+2; 

%  denom2  is  only  used  to  make  sure  the  denom  is  taking  on  a  good  value. 
%denom2  =  zeros (N,N) ; 

PHI  =  zeros (N,N) ;  %  preallocated  for  speed 
for  j  =  1 : N, 


for  k  =  1:N, 

terml  =  (exp(— lj  *  2  *  pi  *  (k— 1)/N)  —  1)  *  sXfft (j,  k) ; 
term2  =  (exp(— lj  *  2  *  pi  *  (j— 1)/N)  —  1)  *  sYfft (j,  k) ; 
denom  =  4  *  (sin  (pi  *  (j  —  1)/N)~2  +  sin  (pi  *  (k—  1)/N)~2); 
if  denom  <  0.000000005, 


%  this  denominator  is  getting  to  be  a  small  value,  but  we 
%  cannot  divide  by  zero.  So  this  checks  to  prevent  that  from 
%  occurring.  Choosing  this  threshold  is  based  on  the  N 
%  dimension.  For  larger  N,  the  threshold  must  be  smaller. 
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%  This  has  major  changes  on  the  outcome  of  the  reconstruction, 
denom  =  0.000000001; 

end 

%denom2(j,  k)  =  denom; 

%  This  PHI  equation  is  from  the  Poyneer  paper. 

PHI(j,  k)  =  (terml  +  term2)  /  denom; 

end 

end 

%  The  Poyneer  paper  specifically  states  the  (1,1)  entry  needs  to  be  0,  but 
%  in  this  filter  implementation,  the  0  is  always  there  anyways. 

%PHI (1, 1 ) =0 ; 

%  the  reconstructed  wavefront.  Should  be  completely  real  but  it  does  have 
%  some  small  imaginary  components  from  roundoff  errors, 
phinew  =  ifft2(PHI); 

phimean  =  phi— mean (mean (phi) ) ; 

%  Error  plot 

figure;  imagesc (abs (phimean— phinew) ) ; 

title (’Error  between  estimate  and  actual');  axis  off; 

%  Original  plot 

figure; imagesc (real (phimean) ) ; 

title  (  ' \Phi  original  data  (mean  subtracted)  ');  axis  off; 

%  Reconstructed  plot 

figure;  imagesc (real (phinew) ) ; 

title ( ' \Phi  estimate  data');  axis  off; 

B.2  CheckAdd.m 

function  [  list,  entry  ]  =  checkadd (  inlist,  inentry  ) 

%[list,  entry]  =  CHECKADD (inlist,  inentry) 

%  Selectively  adds  to  unknowns  list 

%  If  variable  inentry  is  contained  in  the  variable  inlist, 

%  the  output  is  an  updated  entry  with  the  correct  .number 
%  and  the  list  is  not  changed. 
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o, 

o 

%  If  the  variable  inentry  is  not  in  the  variable  inlist,  it  is  added  and  a 
%  new  number  is  assigned.  The  updated  list  and  entry  are  the  output. 

len  =  length  ( inlist ) ; 
add  =  1; 
temp  =  — 1; 

for  j  =  1 : len, 

if  ( inlist ( j ). axis  ==  inentry . axis )  &&  ... 

(inlist  (j).j  ==  inentry. j)  &&  ... 

(inlist  (j).k  ==  inentry. k) 

%  this  entry  already  exists,  so  don't  add  it. 
add  -  0; 

temp  =  j;  %  set  the  output  to  the  existing  entry 

end 

end 

if  add  ==  1, 

%  add  the  unknown  number 
inentry . number  =  len+1; 

%  set  the  list  to  include  the  new  entry 

%  the  .sign  is  removed  as  it  is  not  important  for  the  list  purposes, 
list  =  [inlist,  rmfield ( inentry,  'sign')]; 

else 

%  list  is  not  updated  to  include  new  entry 
list  =  inlist; 
if  temp  >  —1, 

inentry . number  =  temp; 

end 

end 

%  set  entry 
entry  =  inentry; 
end 
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APPENDIX  C: 
Additional  Figures 


43 


Figure  C.l:  SMT  numbered  actuator  topology  for  a  single  segment. 
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Figure  C.2:  SMT  wavefront  sensor  and  actuator  topology  for  a  single  segment. 
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